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A new method for predicting binding affinity in computer-aided 
drug design 
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Jan-Erik Samuelsson 1 

Department of Molecular Biology, Uppsala University, Biomedical Centre, 
Box 590. S-75124 Uppsala and 'Symbicom AB, Glunten, S-75183 Uppsala, 
Sweden 

A new semi-empirical method for calculating free energies 
of binding from molecular dynamics (MD) simulations is 
presented. It is based on standard thermodynamic cycles and 
on a linear approximation of polar and non-polar free energy 
contributions from the corresponding MD averages. The 
method is tested on a set of endothiapepsin inhibitors and 
found to give accurate results both for absolute as well as 
relative free energies. 
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Introduction ■ ,( ; , p\ \ ; ^' 

Computational chemistry ana" molecular modelling have become 
an essential part of the modern dru^ design process. This type 
of methodology is now being used as a complement to experi- 
mental biochemical studies and 3-D structure determination by 
crystallographic or NMR methods. , The most important applica- 
tions of computational modelling in drug design comprise 
(i) methods for finding new Mead compounds', (ii) interactive 
computer graphics for modifying arid manipulating the chemical 
and geometrical structure of inhibitors and (iii) subsequent energy 
and structure refinement using molecular mechanics (MM) lor 
dynamics (MD) calculations [for reviews, see Cohen et al. (1990) 
and Dixon (1992)]. Atpresent,,hpwever, molecular, modelling 
methods can only provide ra&er qualitative information on the 
affinity of different compounds fpr'^a given target site. At best, 
one can propose some tentati ve binding candidates and perhaps 
exclude some that do not seem to match the target site well. This 
may be done by considering molecular shape and electrostatic 
properties, but the quantitative discrimination between different 
ligands in terms of actual binding constants must usually be left 
to the experimental work of organic chemists. The reason for 
this is that it is extremely difficult to* theoretically predict relative 
or absolute binding constants, i.e. free energies, for all but the 
very simplest cases. The only rigorous approach to this type of 
problem today is the so-called free energy perturbation (FEP) 
technique (for reviews, see Beveridge and DiCapua (1989), 
Jorgensen (1989) and Straatsma and McCammon (1992)] which 
is more or less limited to * small perturbations', as will be 
discussed below. 

If we consider a typical case where one wants to determine 
the relative free energy of binding between two compounds* A 

A{w) - B{w) 

AGbwM) I * *<KJLB) 

A(p) - Bip) 
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and B, the problem is described by the thermodynamic cycle 
where AAG^ and AAG^ denote the differences in solvation 
energy between A and B in water and in the (solvated) protein 
site, respectively, and the AG^s are the corresponding binding 
energies; AG^fl) - &GbmM) = AAG&, - The 
same cycle can be used to deterrnine the absolute binding constant 
of B if we let A denote a nil particle. In that case AG^iA) m 
0 and the binding energy is obtained as the difference between 
the absolute solvation energies for B in water and in the protein. 
With the FEP approach one calculates the free energies associated 
with the two unphysical paths A(w) — B(w) and Aip) — Bip), 
corresponding to a mutation of A into B (or the creation of B 
in the case where A is a nil particle). MD (or Monte Carlo) 
simulations are used to collect ensemble averages along the paths, 
which must be rather fine grained in order for the free energy 
to converge properly . If we consider, for example, two enzyme 
inhibitors, the path connecting them will involve changes in the 
molecular charge distribution as well as the creation/anmilation 
of atoms. In [ particular, the,latter jtype of process converges slowly; 
and, thus, if the A ami B i sta& , ; are very different it will be , 
difficult to obtain convergent free energies even with present day 
computer resources (Mitchell and McCammon, 1991; Pearlman 
andKollman, 1991; Straatsma arid McCammon, 1991). It should 
also be noted that most of the computing time in FEP is actually 
spent on uninteresting configurations that correspond to a mixture 
of A and B. Moreover, if large conformational changes are' 
involved this always poses a major problem. Hence, the method, 
only really works' well when # arid B are rather similar and, ui 
{act, all applications related to drug design that have been reported" 
fall within this category (e.g. Wong arid McCammon, 1986; Bash 
etal; 1987; Brooks, \9S9'MeA^KoM^ 
and Brooks, 1991 : Rao and Singh", i$9i ; Tropsha and Hermans, 
1992; the largest 'r^rtuifcauon' to?wJhdch the FEP method has 
been applied seems to be a hydrogen to nexyl group transforma- 
tion reported by Merz 'et al r 199l£ prug design in practice, 
however, often deals with relatively large inhibitors that differ 
considerably from each other so that neither relative nor absolute 
free energies can be obtained with the FEP/MD method within 
reasonable computing time. This dilemma prompted us to try 
to devise a simplified approach that does not involve the 'muta- 
tional paths 1 which is what mainly hampers the FEP method. 

Outline of the method 

As the starring point we take the linear response approximation 
for electrostatic forces, which for polar solutions will yield 
quadratic free energy functions d^otentials of mean force) in 
response to changes in electric fields. This is, for example, the 
familiar result from the Marcus' theory of electron transfer 
reactions (Marcus, 1964). For a system with two states, A and 
B, given by two potential energy functions V A and V B one 
obtains, within the approximation of harmonic free energy 
functions of equal curvature, the relationship (see Appendix and 
for example, Lee et al. (1992) and references therein] 

X = <V S -V A > A -AG AB = <y A -V B > B + AG AB (1) 

385 



J.Aqvfet, C .Medina and J.-E.Sanradsson 




IS 



■ •> if!" 



Fig. 1. Illustration of how the approximation of equation. £2) can be lisetTto estimate the electrostatic contribution to AG lol in a given environment, that can;; 
either be water or protein. The state A corresponds to a system with the/ isolated inhibitor in the T gas phase and a ready made van der Waals cavity in the ; 1 
condensed system. State B is simply the solvated inhibitor, in water or in die protein's active site. These states are given by the two potentials % and V B . Vpr 
any given configuration of the system, equation (2) will reduce to Wi VfL, since the solvent configuration in state A is uncorrelated with the charge ; ' 
distribution of die solute in state B. , . \ ^ , 

calculations on simple systems that corroborate the approxima- 
tion of equation (2). Tests were, for example, carried ou]0y 
comparing the free energy obtained from FEP/MP simulations 
of charging Na + and Ca 2+ ions in a spherical water system 



where AG^j is the free energy difference between B andDl^V 
the corresponding reorganization energy and' <> f - denotes a 
mean evaluated near the niinimum of the potential !. We thus have 



'AG^ = Vfe «AV> A +^AV>i) J ; 



where A V now denotes the energy difference V^y^'V^^/t'; 
consider the hydration of a single ion this'can be stown ^give' 
AGSqj = V4 <VfL s >, i.e. that the electrostatic contribution 
to the solvation energy equals half of the conespoiujing 
ion-solvent interaction energy (Warshel and Russell, 1984;'R6ux 
et a/., 1990). Returning now to our inhibitor binding problem 
we can exploit this result as indicated in Figure t. For each 
solvation process, i.e. solvation of the inhibitor in water and inside 
the protein, we consider two states where the first has the inhibitor 
molecule in vacuo and a non-polar cavity (given, for example, 
by a Lennard-Jones potential) already made in the given environ- 
ment. The second state corresponds to the intact inhibitor 
molecule surrounded by water or the solvated protein. The linear 
response approximation will then again give us AG& = Vi 
< V?- f > where V]l s is the solute -solvent electrostatic term. 
Hence, the electrostatic contribution to the binding free energy 
can be approximated by AG&d = V& < AKfi,> (where the A 
now refers to the difference between protein and water) which 
is obtained from two MD simulations of the solvated inhibitor 
and of the inhibitor— protein complex. 

The validity of the linear response results in the case of ionic 
solvation has been confirmed, for example, in the study by Roux 
et al. (1990). We have also performed some additional 
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S2)V (Aqvist, 1990) with the corresponding ,<V?L' S >, \ft$m>75;$&.. , 
■m : MD trajectories. This yielded fectors relating AG& to . 



of 0.4? for Na + and 0.52 for Ca 2+ ybpUv values. r^m^cJo^tp, . . 
me predicted result of Vi. A siir&aVtest^ 
a methanol molecule, given by the OPUS potential (Jorgeijseii^ 
1986), in water gave a AG&: < ratio of 0.43. Some 
calculations were also done on rigid two-atom dipoles with the 
charge centres separated by 4 A. For a Q = ±2 (elementary 
charge units) dipole a ratio of 0.53 was obtained while Q = ± 1 
and Q = ±0.25 gave ratios of 0.47. and 0.49 respectively. 
Charging a sodium ion given by the parameters in Aqvist ( 1990) 
to +0.25 charge units gave a ratio of AG&: < VfL s > = 0.55. 
Hence, the approximation of equation (2) seems to be reasonable 
although it might be somewhat less accurate tor very weak fields. 
Whether this is due to non-linear effects or simply to the fact 
that the relative fluctuations cf Vf_ s are larger is not entirely 
clear and deserves further study. 

Compared to the electrostatic problem above it is a more 
difficult question how to account for the contribution of non-polar 
interactions and hydrophobic effects to the free energy of binding, 
which we will call AGjfe. What is desirable is to be able to 
estimate this contribution from the non-polar (or van der Waals) 
interaction energies. The liquid theories of Chandler and 
coworkers (Pratt and Chandler, 1977; Chandler et al s 1983) have 
been successfully used to analyse hydrophobic effects and to 




2 4 6 8 10 

? . ■ ■ <■ rr: . 



length of carbon chain , : ; % 

Fig. 1. Calculated dependence of the mean solute -solvent interaction 
energy on the length of the carbon chain for n-alkanes solvated in water. 

calculate free energies of transfer for some non-polar molecules 
(Pratt and Chandler, 1977), but no analytical treatment of that 
kind seems possible for solvation in an inhomogenous environ- 
ment such as a protein's active site. It can, however, be noted 
that the experimental free energy of solvation for various 
hydrocarbon compounds, such as r^alkanes, depends 
approximately linearly oh the length of the carbon chain both 
in their' own liquids as well as in water (Ben-Nairn and Marcus, 
,1984). We have canied.out^MD siniulationsof /f-alkanes solvated 
in Vater (Figure 2) *imd;ura non-polar van der Waals solvent 
(not shown) which mdicate that (he mean solute -solvent 
t interaction energies alko^vary approximately linearly with the 
number of carbons in chain (the relationships being different 
in different solvents;, of course). It would thus seem possible 
that a simple linear approximation of AGJ^ from <AKJ?J > 
might be able to account for the non-polar binding contribution. 
For instance, if we let a be some appropriate measure of the size 
of the solute and if the solute - solvent van der Waals interaction 
energies and the corresponding non-polar free energy contribu- 
tions (both in water and protein) depend linearly on o\ such that, 
< Vf»> = oyy, <V*"> - ajj. AGf" = /S^aand AG^ 
= ft/j, then we obtain AGJ& = &p ' &w <A^ >. Since 

it seems difficult to derive a factor relating the two quantities 
in a reliable way from purely theoretical consideratioas. we take 
the approach here to empirically try to determine such a 
relationship that is capable of reproducing experimental binding 
data. Thus, we will approximate the free energy of binding by 

AGferf = Vi <Wf. s > + oc <A^> (3) 

and attempt to determine the parameter a by empirical calibration. 

As a test system we have chosen endothiapepsin (EP) belonging 
to the family of aspartic proteinases (see, e.g. Fruton, 1976; 




Fig. X piemical. strucmres/of the inhibitors /,,.... / 3 used in the 
c^aaations: s y ,' f f ^/v>Y# * ' ' V 

Davies, 1990), a class of enzymes for which numerous studies 
of inhibitor binding have -been reported. Crystal structures of 
native endothiapepsin and inhibitor complexes have been 
published by Blundell and coworkers (Foundling et al, 1987; 
Veerapandian et al., 1990). In this work we use the crystal 
structure of endothiapepsin in complex with the inhibitor H2 18/54 
(Figure 3, / ( ), recently determined by Symbicom AB (D.Ogg 
et aL, in preparation), as the structural starting point for our 
computations. Four other inhibitor compounds were used in this 
work and their chemical structures are also shown in Figure 3. 
For simplicity we label these molecules i } / 5 . Experi- 
mental binding data for the inhibitors have been obtained by 
Falknas and Deinum (A. Hassle, unpublished results). 

Computational details 

Starting from the experimental structure of the EP-/ t complex 
the other inhibitors were manually built into the EP active site 
using the FRODO (Jones, 1978) and Hydra (Hubbard, 1984) 
graphics software. This model building appears relatively 
straightforward since all the compounds contain the invariant 
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Fig. 4. Stereo view of the 50 ps mean MD structure of the EP-/ ( complex (thin 
etal., m preparation). The active site of the enzyme is shown with the bound ' 

hydroxyethylene transition state isostere adjacent to a peptide 
group, which was also present in the inhibitors studied by 
Blunders group. After initial energy minimization, MD 
simulations consisting of 25 ps of equilibration and 50 ps of data 
collection were performed for each of the inhibitors, solvated 
in water as well as bound to the solvated protein. 

A modified version (J.Aqvist, unpublished) of the ENZYMDC 
program (Warshei and Creighton, 1989) was used for all MD 
simulations together with the GROMOS potential (van Gunsteren 
and Berendsen, 1987). This force field was revised by us with 
respect to the oxygen-extended carbon (CH„) interactions, since 
the standard GROMOS parameters were found from PEP simula- 
tions to give an incorrect .'value of AG^ = +0.1 ± 
0.5 kcal/mol for methane (J.Aqvist, unpublished results). By 
changing the oxygen repulsive Lennard- Jones 6/12 parameter 
(4 0 ) for O-CH, interactions from 421.0 to 793.3 (kcal/mol) 
A 6 , AGhydr for methane becomes 2,5 ± 0.4 kcal/mol, in 
better agreement with the experimental value of 2.0 kcal/mol 
(Ben-Nairn and Marcus, 1984). Note that this revision is quite 
important for a correct dewription of the hydrophobic effect. Ako 
the Lennard -Jones parameters for charged carboxylate groups 
■have been revised earlier in order to better reproduce experi- 
mental solvation energies (Aqvist etal, W^V'The present 
calculations used a spherical water droplet of radius 16 A for 
'the simulations of inhibitors in solution and a corresponding (16 
A sphere containing both protein atoms and water in the simula- 
tions of bound inhibitors. The atoms within this sphere were free 
to move while protein atoms beyond 16 A were restrained to their 
crystallographic positions. An mteracuon cut-off ? radius of § A 
was used, except for interactions involving the inhibitor and the 
MD time step was 0.001 ps. The equilibration phase of the 
protein simulations consisted of 5 ps of successive heating of the 
system and weakening of harmonic positional restraints that were 
applied to the protein atoms. After this period all restraints within 
the 16 A sphere were set to zero and the system was further 
equilibrated at the final temperature of 298 K for another 
20 ps. The r.m.s. coordinate deviation for the inhibitor atoms 
between the following 50 ps average MD structure and the 
experimental EP— /f complex is 0.94 A. Figure 4 shows a 
superposition of these two structures and it can be seen that the 
agreement is quite ; 



lines) superimposed on the corresponding crystallographic structure (D.Ogg 
inhibitor in the centre of the picture. 



Results and discussion 

We used the first four inhibitors (/j, . . ., / 4 ) as our calibration 
set in order to determine the optimal value of the parameter a 
and examine the success of equation (3) in reproducing the 
experimental results. Table I shows the observed and calculated 
absolute free energies of binding for the inhibitors I h . . I 4 
using the value of a = 0.161 which optimizes the r.m.s. 
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Table I. Calculated and observed (Palknas and Deinum, unpublished results) 
absolute binding free energies for inhibitors f u . . „ / 4 to endothiapepsin 
(in kcal/mol) 






1 

2 
3 
4 


-11.16 -10.69 
-8.23 -7.93 

-11.14 -11.67 
-6.29 -6.57 




Table u\ Calculated and observed, relative binding free energies for 
inhibitors / ( , . . .', / 4 to endothiapepsin (in kcal/mol) ^ 







1-2 


-2.92, 


: -2.76 


1-3 


-o.oi' 


+0.98 


1-4 


-4.87 


-4,12 


2-3 


+2.91 


+3.74 


2-4 


-1.95 


-1.36 


3-4 


-4:»6 


-5.10 . 



agreement with the experiment;' This value of « gives a mean 
unsigned error of 0.39 kcal/mol for the calibration set and the 
largest error- is 0.53 kcal/mol .for inhibitor fj. 1 Initially^we had 
not expected such an accurate tWor-the absolute binding' energies 
but had mainly hoped to be able to apply and calibrate equation 
(3) with respect to the relative; ones. If me calibration instead 
is done using relative binding energies one obtains a-very similar 
value of a (0. 169) with mean unsigned errors that are virtually 
identical to those above. This result is quite remarkable and, thus, 
lends further support to the approximation of equation (3). The 
calculated and observed relative free energies of binding (using 
<* = 0.161) are shown in Table II and the mean unsigned error 
is 0.59 kcal/mol. The largest error here is 0.99 kcal/mol for the 
/ ( // 3 selectivity, but all other relative free energies have the 
correct sign and are within 0.8 kcal/mol of the experimental 
result. It is particularly interesting to note that the calculations 
are able to discriminate the low-affinity inhibitor / 4 quite well 
from the high-affinity ones (/, and / 3 ). 

Encouraged by the above results, we wanted to try to assess 
the predictive power of our approach uy modelling an inhibitor 
not present in the calibration set, For this we chose the inhibitor 
/ 5 which, as can be seen from Figure 2, differs significantly in 
its chemical structure from any member of the calibration set. 
This molecule was built into the EP active site and subjected to 
the same simulation procedure as the other inhibitors. The 
predicted absolute free energy of binding for J$ is 
-9.70 kcal/mol which is in excellent agreement with the 
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Table UI. Calculated and observed relative binding free energies involving 
inhibitor U (in kcai'mol) 



i-5 


-1.46 


-0.85 


2-5 


+ 1.47 


+ 1.91 


3-5 


- 1 .44 


-1.83 


4-5 


+ 3.41 


+ 3.27 



corresponding observed result of AG^^s) = -9.84 kcal/mol. 
The calculated relative binding free energies* with respect to the 
four inhibitors in the calibration set are given in Table HI, where 
it can be seen that all the pairwise selectivities involving / 5 are 
correctly predicted by our simulations, the maximum error in 
this case being 0.61 kcal/mol. The above results indicate that 
the simple linear approximation of the free energy given by equa- 
tion (3) is able to describe the main physics (electrostatic and 
hydrophobic interactions) of the binding process in a satisfactory 
way. (In retrospect, it is also of interest to see what the effect 
would be of choosing different inhibitor subsets for the calibra- 
tion of a. It then turns out that the possible subsets of four 
inhibitors give values in the range 0.158 < or < 0.165. For 
example, a = 0.165 is obtained if/, is left out of the calibration 
and its predicted binding energy then becomes - 1 1 .31 kcal/mol. 
If all five inhibitors are included in the calibration one obtains 
a value of ce = 0.162.) *' '■ " \7 

It is, of course, important to try to identify and separate the 
different ; types of errors involved in the present method, and 
suggest how they possibly could be dealt with. There are basically 
four sources of errors, namely (i) inaccurate starting structures 
due to incorrect model building, (ii) possible deficiencies in 
potential energy functions, (iii) poor MD convergence, for 
example, due to short trajectories and (iv) the approximation of 
equation (3) itself. The only remedy for errors of the first type 
[ is to obtain as much structural information as possible and in 
order to assess their magnitude it would be desirable to, cany 
out calculations on a set of inhibitor complexes whose 3-D 
^tmcmres^haye all been experimentally determined. Although 
^MM/Mp^rwtentials; are continuously being refined by .many 
' rese^h' groups errors of the second type may still be consider- 
,.able.Jt;is,' f^rixample, not yet entirely clear how well customary 
protein *fprce fields actually reproduce relevant . energetic 
' prpperties; In'our case we can note that the force field used here 
has required some revision in order to reproduce essential 
solvation properties (see above). 

The convergence properties of MD simulations depend on how 
far from equilibrium the initial structure is, but judging from the 
present study it seems that one can reach satisfactory convergence 
within reasonable computing time. For example, by comparing 
means over the first and second halves of the MD trajectories 
we obtain mean (over all five inhibitors) errors of ±0.35 and 
±0.75 kcal/mol for V?? s and vf. s respectively, in the 
protein; the corresponding errors in water are ±0.46 and 
±0.62 kcal/mol. This would yield a nominal error rar. 0 e of 
±0.82 kcal/mol in equation (3) originating from the MD 
convergence uncertainty. Here, it can be noted that an advantage 
with the present approach is probably that it focuses mi simulation 
of the thermodynamically relevant states compared, for example, 
to FEP calculations where most of the computing time is spent 
on the paths between such states. The errors associated with the 
approximation of equation (3) are at the present stage difficult 
to estimate quantitatively although the agreement with experi- 
mental binding data obtained here indicates that the linear 



approximation is reasonable. A larger calibration set would 
obviously be desirable but this may be regarded as a matter of 
refinement, which we leave for future work. As mentioned above, 
Ben-Nairn and Marcus (1984) have shown for several classes 
of hydrocarbon chain containing organic compounds that a linear 
fit o^ A<7 W / = kn + /, where n is the number of carbon atoms 
in the chain, is quite accurate both in water and non-polar 
solvents. It can, however, be noted that in some cases ihe 
extrapolation, /, of AG^ to zero chain length is non-zero 
(e.g. for n-alkanes in water / = 1.42 kcal/mol). This might 
suggest that an additional constant term should be added to 
equation (3), reflecting a difference in extrapolation to 'zero 
inhibitor size' (between the water and protein environments) of 
the non-polar contribution to AG^. However* such a term 
would again be difficult to derive in a reliable way for the systems 
we are dealing with and it does not seem warranted at the present 
stage to introduce another empirical parameter in equation (3) 
(in fact, doing so does not improve the fit reported in Table I 
significantly). H is also important to emphasize here that equation 
(3) with the above parametrization of a is not an equation for 
the individual solvation energy terms, since the factor a represents 
the combined effect of several energy/free energy relations 
[however, one could of course try to parametrize equation (3) 
for solvation energies which* is of interest for calculating 
free energies of transfer between different solvents]. The fact 
that we always deal with solvation energy, differences fhay;also 
cause some cancellation of possible systematic errors. Further^ 
more^due to the larger weight given by equation (3) ; to the 
&Vfj s term one might perhaps expect the formula 'to give 
better results for polar inhibitors, as long as the electrostatic linear 
response approximation holds. When comparing inhibitors of 
different ^charge states long-range corrections of the Born type 
(see, e.g. Srraatsma and Berendsen, 1988; Aqvist, 1990) would 
obviously become important and may make accurate predictions 
more difficult, although counter ions may be used to keep the 
system overall neutral. \ - 
It also remains to be seen whether the empirical parameter ot 
is readily transferable to other systems or whether it will display 
.some system-dependency. One, thing that is clear, given the fact 
■that the : panntietrization of force fields can differ considerably f 
. is that ajs.likeiy to be force field-specific. If we are unfortunate 
s this parameter might also, for example, reflect the nature, of Uie - 
particular; protein site. But, even if this turns out to be the case 
we* feel that the present approach could be^uite efficient for 
discriminating between a large number of potential inhibitors to 
a given target site, thereby reducing the need to synthesize all 
of them. One would then, of course, have to go through an initial 
round of calibration but this might be worth the effort. It should 
also be emphasized that there does not really exist any reliable 
computational method today that could match the qualitative 
level of results reached here for such large and widely differing 
inhibitors. 
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Appendix 

Figure 5 depicts schematically the free energy functions for a 
system with two states defined by potential surfaces V A and V B . 
The free energy functions are defined as the potential of mean 
Torte with respect to some reaction coordinate, X, that is, the 
probability distribution for X on the potential surface V A is 

P „(X) = <tur-x\> A = i 8( y~Sr^' dr < 4 > 

and the corresponding free energy function is then defined as 
A&CX) = -r'hipi(» (5) 

where 0"' = k B T, Analogously, the free energy function on the 
V B surface is written as 

&8dX> = -J3*' InprfX) + AG AB (6) 




Fig. 5. Schematic free energy curves for a system obeying the linear 
response approximation. Sec text for details. 

where the constant AG AB is the (reaction) free energy associated 
with changing the state of the system fromv4 to B. In the present 
case, we consider this change to only involve electrostatic 
variables so that the difference between V A and V B reduces to 
coulombic terms. The reaction coordinate, X, describes the 
change from state A to state B and can be chosen in different 
ways (see, e ; g. King and Warshel, 1990), for instance, as the 
energy gap between the two surfaces or as- a suitable charging 
parameter . The essence of the linear response approximation is 
that the free energy functions assume parabolic form and are 
characterized by the same force constant (i.e. the system responds 
linearly to electrostatic forces with a certain dielectric constant). 
This is the standard dielectric continuum result of which Marcus' 
ET theory provides a nice illustration (Marcus, 1964). Equation 
(1), which relates the. reorganization energy (X) to the reaction 
free energy and mean .energy gap, then follows from the simple 
geometric relationship between the two free energy curves.; It 
should thus be noted 'that the factor of Vx in equation (2), that 
can be derived in different ways (e.g. Warshel and Russell, 1984; 
Roux et ai , 1990), 'originates from the assumption of. harmonic 
free energy curves ofrequal curvature: ;If the curvatures of the 
parabolas are not equal one obtains an additional term of - ViAX 
in equation (2)^ where the difference in, Reorganization energy 
between the two states is related to the. origin shift and force 
constants of the parabolas by AX - {Ak°Y(k B - k A ). 

An alternative derivation of equation (2) can be carried out 
starting from the perturbation expression for the free energy 
difference, which can then be written as (Zwanzig. 1954) 

AG = -jtr 1 in <c"^« ~ V *'> A (7) 
Expansion of the exponent and logarithm then yields 

AG = -jrMn<l-0^-K<) + £) 

= -j3-'ln|l-j3 <[V B -V A )> A + 4) 

«Vb-Va?>a 

- fc<-$ < y B -y A >4 + 4) 

<(^-K^>J ; + ■ I (8) 
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which can be rearranged as 

AG * <A(^>^ - (|) <{AV-<AV> A ) 2 > A + 
<£) <(AK-<AK> H ) 3 > y4 

6 

- (£)\<(AV-<AV> A )*> A - 

24 

3<(AK-<A^> /1 ) 2 >/1 + . . . (9) 

where AK = K^-l^. In the same way we obtain by averaging 
on the potential surface V B 

AC = <AV> B + (|) <(AV-<AV> B f> B + . . . (10) 

Adding the two equations then gives us 

AC = Vi\<AV> A + <AV> B \ 

- (*) [<(AV-<AV> A ) 2 > A - 

<{AV-<AV> B f> B \ + . . (11) 

Although equation (2) could now be retained by just dropping 
the terms of second and higher order, it should be pointed put 
that a simple truncation to first order in AV would correspond 
to the approximations AG < AV> A (from equation 9) and 
AG f <AV> B (from equation 10), that is, an assumption of 
theanean energy gaps being equal on the two surfaces, which 
is clearly not justified. Therefore, the form of equation (2) is \ 
formally obtained by truncating terms of third order and r higrjer 
and by assuming that the mean square fluctuations of the energy 
gap on the two surfaces are equal and, thus, will cancel. Hence; 
the approximations here are identical to those behind Figure 5, 
that is, truncation of the energy gap expansions beyond the second 
moment is equivalent to the assumption of harmonic free energy \- •* 

functions and the cancellation of the fluctuation terms is equivalent ■ k * 
to the assumption of equal force constants for the two states (the ' • . 

proof of this is straightforward). In order for equation^) to Be 
.valid it As, actually, sufficient with the somewhat . weaker ^..jl.y. 
. assumption of 1 symmetric free energy functions* with equal; u*v.ttf; 
expansion coefficients. It can also be noted that equations (9) and 
(10) provide a means of obtaining reorganization energies from , 
simulations and, correspondingly, in cases where the equal ' ; , 
curvature approximation does not hold equation (1 1) provides t ;* 
a means for evaluating the difference between the force constants . ' • - 
and, thus, for the difference in reorganization energy. 
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